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In this paper we introduce a model which provides a new approach to 
the phenomenon of stochastic resonance. It is based on the study of the 
properties of the stationary distribution of the underlying stochastic process. 
We derive the formula for the spectral power amplification coefficient, study 
its asymptotic properties and dependence on parameters. 

Introduction 

The notion of Stochastic Resonance appeared about twenty years ago in the 
works of Benzi et al. [l j and Nicolis |2. in the context of an attempt to 
explain the phenomenon of ice ages. The modern methods of acquiring and 
interpreting climate records indicate at least seven major climate changes in 
the last 700,000 years. These changes occurred with the periodicity of about 
100, 000 years and are characterized by a substantial variation of the average 
Earth's temperature of about 107^. 

The effect can be explained with the help of a simple energy balance model 
(for an extended review on the subject see [B]). The Earth is considered as a 
point in space, and its temporally and spatially averaged temperature X(t) 
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satisfies the equation 

2 irt 

X(t) = -U , (X(t))-Qmn{—), (1) 

where U(X) is a double-well potential with minima at 278. 6K and 288. 6K 
and saddle point at 283. 3-K" and wells of equal depth. The second term in (JTJ) 
corresponds to a small variation of the solar constant of about 0.1% with a 
period of T = 100, 000 years due to the periodic change of the eccentricity of 
the Earth's orbit caused by Jupiter. The influence of this term reflects itself 
in small periodic changes of the depths of the potential wells. In this setting, 
the left well is deeper during the time intervals (kT, (k + |)T), whereas the 
right one is deeper during the intervals ((A; + |)T, (k + 1)T), k — 0, 1, 2, . . . 

The trajectories of the deterministic equation (0) have two metastable 
states given by the minima of the wells. Due to the smallness of the solar 
constant Q no transition between these states is possible. In order to obtain 
such transitions Benzi et al. [T] and Nicolis j2] suggested to add noise to the 
system which results in considering the stochastic differential equation 

X^(t) = -U>{X*> T {t)) - Q sin (^) + V~eW t , (2) 

e > 0, W a white noise. 

Now one can observe the following effect. Fix all parameters of the system 
except e and consider the typical behaviour of the solutions of (j2j) for different 
values of e. If the noise intensity is very small, the trajectory only occasionally 
can escape from the minimum of the well in which it is staying, and one can 
hardly detect any periodicity in this motion. If the intensity is very large, the 
trajectory jumps rapidly but randomly between the two wells and therefore 
also lacks periodicity properties. An interesting effect appears when the 
noise level takes a certain value e$: the trajectory always tends to be near 
the minimum of the deepest well and consequently follows the deterministic 
periodic jump function which describes the location of the deepest well's 
minimum. It is very important to note that to produce this effect one needs 
all three of the following components to be present in the system (J2J): the 
double-well potential for bi-stability, the noise to pass the potential barriers, 
and a small periodic perturbation to change the wells' depths. 

The following are natural questions arising in the context of these qual- 
itative considerations: how can one measure periodicity of the trajectories 
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and, consequently, how does the quality of tuning of the noisy output to the 
periodic input be improved by adjusting the noise intensity e? 

The formulation of the latter question suggests to consider the system (J2J) 
as a random amplifier. The random system receives the harmonic signal of 
small amplitude Q and usually large period T as input. The stochastic 
process X e,T (t) is observed as the output. The input signal carries power 
Q 2 at frequency 1/T. The random output has continuous spectrum and 
thus carries power at all frequencies. Benzi et al. pQ considered the power 
spectrum of the output for different values of e and discovered a sharp peak 
at the input frequency for a certain optimal value of Eq. This means that the 
random process X So,T (t) has a big component of frequency 1/T. The effect of 
amplification of the power carried by the harmonic considered as a response 
of the nonlinear system (J2J) to optimally chosen noise was called stochastic 
resonance. 

In the past twenty years more than three hundred papers on this sub- 
ject were published. An extensive description of the phenomenon from the 
physical point of view can be found in jlj and [5]. The notion stochastic 
resonance is now used in a much broader sense. It describes a wide class 
of effects with the common underlying property: the presence of noise in- 
duces a qualitatively new behaviour of the system and improves some of its 
characteristics. 

Although stochastic resonance was observed and studied in many physical 
systems, only few mathematically rigorous results are known. The approach 
of M. Freidlin is briefly outlined in the next section of this paper. In sec- 
tions 121 El and 0] we introduce discrete-time Markov chains with transition 
probabilities chosen in such a way, that on a large temporal scale the attrac- 
tor hopping behaviour of the underlying diffusion process is imitated in the 
limit e —>■ 0. We investigate stochastic resonance for the Markov chains. The 
last section is devoted to generalizations and discussion. 

1 Large deviations approach 

In this section we briefly survey rigorous mathematical results obtained by 
M. Freidlin in jH] using the theory of large deviations for randomly perturbed 
dynamical systems, developed in Freidlin and Wentzell (see (7j). Though the 
results of are valid in a quite general framework, we confine our attention 
to a simple example of a diffusion with weak noise. 
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Consider the SDE in R 

X^{t) = -U\X^{t\ t -) + ^eW{t\ (3) 

where W is a white noise and U'(x, t) = J^U (x, t), with a time dependent po- 
tential just periodically switching between two symmetric double well states, 
i.e. 

U(x 1 t) = ^U(x)l [kM i ) (t) + U(-x)l l ^i M1) (t) 1 

where U(x) has local minima in x = ±1 and a saddle point in x = 0, 
lim^i^oo U(x) = oo. We also fix the depths of the wells by two numbers 
< v < V, assuming that U{-1) = -V/2, U{1) = -v/2, and U(0) = 0. 
Note, that X e,T is a Markov process which is not time homogeneous. In 
the following Theorem time scales are determined in which some form of 
periodicity is observed. 

Theorem 1 Suppose T = T{e) is given such that 

limelnT(e) = A > 0. 

a) If A < v, then the Lebesgue measure of the set 

{t E [0, 1] : \X e ' T{£ \T(e)t) - sgnX | > 5} 

converges to in Px probability as e — > 0, for any 5 > 0. 

b) If A > v, then the Lebesgue measure of the set 

{tE[0,l] : |X £ ^)(T( £ )t)-0(t)|>5} 
converges to in Px probability as e — > 0, for any 5 > 0, where 

<f>(t)=^-l lkM i ) (t) + l l ^i M1) (t) 

k>0 

and Px denotes the law of the diffusion starting in X . □ 

It is nessesary to explain why A = v is critical for the long time behaviour of 
the diffusion. At least intuitively, the answer follows from the asymptotics of 
the mean exit time from a potential well for the time-homogeneous diffusion. 
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If the diffusion starts in the potential well with the depth v/2, its mean time 
E(r(e)) needed to leave the well satisfies 



elnE(r(e)) -»• v, e^O. 

according to Freidlin and Wentzell [Zj. This means, X £,T ^ can leave nei- 
ther the deep well with the depth V/2 nor the shallow one with the depth 
v/2 in time T(e) of order e A / £ if A < v. Therefore, X £ ' T ^ stays in the 5- 
neighbourhood of the minimum of the initial well. On the other hand, if 
A > v , X e,T ^ has always enough time to reach the deepest well. In both 
cases, the Lebesgue measure of excursions leaving the 5-tube of the deter- 
ministic periodic function is exponentially negligible on the time scale T(e) 
as £ — > 0. 

The Theorem suggests the time scale which induces periodic and deter- 
ministic behaviour of the system (JHJ), and the Lebesgue measure as a measure 
of quality. In fact, it only gives a lower bound for the scale. In the next sec- 
tion, in the framework of discrete Markov chains approximating the diffusion 
processes just considered, we investigate different measures of quality which 
provide unique optimal tuning. 

2 Markov chains with time-periodic transi- 
tion probabilities 

For m G N, consider a Markov chain X m = (X m (k))k>o on the state space 
S = {—1, 1}. Let P m (k) be the matrix of one-step transition probabilities at 
time k. If we denote vr~(/c) = P(X m (k) = —1), vr+(/c) = P(X m (k) = 1), and 
write P* for the transposed matrix, we have 

*m( k + 1 ) \ _ p*(l,\ ( *m( k ) 

In order to model the periodic switching of the double-well potential in 
our Markov chains, we define the transition matrix P m to be periodic in time 
with half-period m. More precisely, 



Pm(k) 



Pi, < /c(mod 2m) < m — 1, 
P2, m < fc(mod 2m) < 2m — 1, 
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with 

1 V V 

1 — px px 



Pi 




qx v l-qx l , , 
1 — qx v qx v s 
px v 1 — px v 



where <p = (p(e, p, V) = pe~ v ^ £ , ip — ip(e, q, v) = qe~ v ^ £ , x = e _1 / £ , < p, q < 
1, < v < V < +oo, < e < +oo. Sometimes, it will be convenient to 
consider x G [0, 1]. In these cases the ends of the interval will correspond to 
the limits e — > and e — > oo. 

In this setting, the numbers V/2 and v/2 clearly have to be associated 
with the depths of the potential wells, e with the level of noise. According to 
the Freidlin-Wentzel theory, the exponential factors in the one-step transition 
probabilities just correspond to the inverses of the expected transition times 
between the respective wells for the diffusion considered in the preceding 
section. This is what should be expected for a Markov chain in equilibrium, 
modulo the phenomenological pre-factors p and q. They model the pre- 
factors appearing in large deviation statements, and add asymmetry to the 
picture. 

It is well known that for a time- homogeneous Markov chain on S with 
transition matrix P one can talk about equilibrium, given by the stationary 
distribution, to which the law of the chain converges exponentially fast. The 
stationary distribution can be found by solving the matrix equation n = P*n 
with normalizing condition 7r~ + 7r + = 1. 

For non time homogeneous Markov chains with time periodic transition 
matrix, the situation is quite similar. Enlarging the state space S to S m = 
{— 1, 1} x {0, 1, ... , 2m— 1}, we recover a time homogeneous chain by setting 

Yjk) = (X m (k), k{mod 2m)), k > 0, 

to which the previous remarks apply. For convenience of notation, we assume 
S m to be ordered in the following way: 

S m = [(-1, 0), (1, 0), (-1, 1), (1, 1), ... , (-1, 2m - 1), (1, 2m - 1)]. Writing 
A m for the matrix of one-step transition probabilities of Y m , the stationary 
distribution Q = (q(i,j))* is obtained as a normalized solution of the matrix 
equation (A^ — E)Q = 0, E being the unit matrix. We shall be dealing with 
the following variant of stationary measure, which is not normalized in time. 
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Definition 1 Let 7r m (k) = (n m (k), n+(k))* = 2m(g(-l, k), q(l, k))*, < 
k < 2m— 1. We call the set7c m = (7T m (A;))o<fc<2m-i the stationary distribution 
of the Markov chain X m . 

The matrix A m of one-step transition probabilities of Y m is explicitly given 
by 



A, 
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A m has block structure. In this notation means a 2 x 2-matrix with all 
entries equal to zero, Pi, and P2 are the 2- dimensional matrices defined in (£Q). 

Applying some algebra we see that (A^ — E)Q = is equivalent to 
A' Q = 0, where 



/ P-E 
P* —E 



A„ 













V 





















-E 



p* 





—E 

p* 
r 2 



\ 






- E J 



and P = P^P; 



P* = (P*) m (P*) m . But A' m is a block-wise lower diagonal 
matrix, and so A' m Q = can be solved in the usual way to give 

Theorem 2 For every m > 1, the stationary distribution ir m of X m with 
matrices of one-step probabilities defined in ^ is: 



+ 



(p + 4) (p + ^ 1 + (1 - (p - ipy 
if (p — ip (1 — tp — ipy 



ir + (l) 
ir~(l + m) = 7T+(/), 

7r+(Z + m) =tt-(Z), 0</<m-l. 



(5) 
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Proof: vr m (0) satisfies the matrix equation ((P 2 *) m (P*) m -P)vr m (0) = with 
additional condition vr m (0) + vr+(0) = 1. To calculate (P 2 *) m (P 1 * 
formula for the m-th power of 2 x 2-matrices, which results in 



/ 2 -p_i _i -pi i V 1 l-P-i-i 

(p-i,-i + pi,i - ir / i - p-i,-i -(i - p-l-i; 



+ 2-p_i_i-pi jl V-( 1 ~Pi,i) 
Using some more elementary algebra we find 



(P 2 *) m (P*) m = (pmp ; 



m\ * 



1 — if) if) \ / 1 — ip (f 
(p 1 — p J V ^ 1 — if) 



■v 3 + ^ V V' / y? + ?/> \ 1 - 1 - 



v? + ^ V ^ 

from which a straightforward calculation yields 

(p + if){i — cp — ipy 



m( } (^ + ^)(1 + (1-V5-V) m )' 

+rm = y> + - - y>r 

7Fml j (v? + ^)(l + (l-v?-^) m )' 

To compute the remaining entries, we use ir m (l) = (P 1 *) / 7r m (0) for < I < 
m-l, and ir m (l) = (P 2 *)'(P 1 *) m vr m (0) for m < I < 2m - 1 to obtain ©. Note 
also the symmetry 7r~(/ + m) = vr+(/) and 7r+(/ + m) = 7r~(Z), < / < m — 1. 
□ 



3 Spectral power amplification 

The chain X m can be interpreted as amplifier of a signal. Our stochastic 
system may be seen to receive a deterministic periodic input signal which 
switches the double depths of the potential wells in (jl]), i.e. 



I'm (/) 



V, < /(mod 2m) < m — 1, 
v, m < /(mod 2m) < 2m — 1. 

The output is a random process X m (k). 
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The input signal I m admits a spectral representation 

^ 2m- 1 



2m-l 

. , 2nik , 



, c n ,((i u 2 ™ a 

2m 

a=0 



where c m (a) = (l/2m) ^S)" 1 hi^e^' 1 is the Fourier coefficient of frequency 
a/2m. The quantity |c 2m (a)| 2 measures the power carried by this Fourier 
component. We are only interested in the component of the input frequency 
l/2m. Its power is given by 

In the stationary regime, i.e. if the law of X m is given by the measure 7r m , 
the power carried by the output at frequency a /2m is a random variable 

1 2m- 1 

<-W = 5S £ *-(')««'■ 

2=0 

We define the spectral power amplification as the relative expected power 
carried by the component of the output with frequency 

Definition 2 The spectral power amplification coefficient of the Markov 
chain X m with half period m > 1 is given by 



Vr, 



|e*»(6»(i))I s 



|c m (l)| 2 ' 

Here E nm denotes expectation w.r.t. the stationary distribution n m . 

The explicit description of the invariant measure now readily yields the 
following formula for the spectral power amplification. 

Theorem 3 Let m > 1 . The spectral power amplification coefficient of the 
Markov chain X m with one-step transition probabilities (J1J equals 

4 (y?-^) 2 



(V-vf (^ + V,)2 + 4(l-^-^sin 2 (^-)- 
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Proof: Using (j5J) one immediately gets 

1 2m— 1 1 n i m—1 

fc=0 k=0 

2 (p-ip ( 1 1 
m^ + V^l-e™ l-(l-v9-V)e- 

Some algebra and an appeal to (JHJ) finish the proof. □ 
Recall now that the one-step probabilities Pi and P2 depend on the pa- 
rameters < p, q < 1 and, what is especially important, on < e < 00 
which is interpreted as noise level. Our next goal is to tune the parameter 
£ to a value which maximizes the amplification coefficient t] m = f] m (e) as a 
function of e. 



4 Extrema and zeros of rj m (e). 

In this section we study some features of the function r) m (e) and its depen- 
dence on m G N, Q < v <V < 00 and the pre-factors < p < 1, < g < 1. 

After substituting e~ l ' £ = x and writing r) m (e) = T} m (x), this function 
takes the form 

/ \ = __4 (PX V -qx v f 

Vm[X> (V - v) 2 (px v + qx v ) 2 + A{l-px v -qx v ) sin 2 (£) U 

In what follows, we assume x G [0,1]. The boundaries x = and x — 1 
correspond to the limiting cases e = and e = 00. Denote a m = esc (^) 2 > 
1, m > 1. 

Our main result on optimal tuning is contained in the following theorem. 



Theorem 4 a) We have rj m (x) > 0, rj m (0) = 0. 

b) Let < (3 = y < 1 and m > 1 be fixed. There exists a continuous 
function 

b(q; m, (3) - a/%; m, (3) 2 - 4a(g; m, f3)(2 - q)q 

P-{Q) =P-{T,P,rn) = — , 

2a{q; m, p) 

where a(q; m, (5) — 1 — a m q(l — 0), b(q; m, (3) = 2 — 3(1 — (3)q + a m (l — (3)q 2 , 
with following properties: 
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i) P-(q) > 0, q G [0, 1] and p_(q)=0<^q = 0; 

ii) p-(q) < q, q £ [0, 1] and P-(q) = q q = or q = 1, m = 1; 



.. dp-(q;m,/3) 
m J 



dq 



P. 



g=0 



Moreover for m > 2 

1) If (p,q) E Uq = {(p,q) : < q < 1,0 < p < P-{q)}, r] m {x) is strictly 
increasing on [0, 1]. 

2) If(p,q) G U\ — {(p,q) : < q < l,p_(q) < p < q}, r] m {x) has a unique 
local maximum on [0,1]. 

3) If (p,q) G U% — {(p,q) '■ < q < 1, q < p < 1}, r) m (x) has a unique 
local maximum on [0, 1] and a unique root on (0, 1]. (See Fig. ) 

c)For any 5 > there exists Mq = Mo(p, q,/3, 5) such that for m > Mq 
the coordinate of the local maximum x m G [x m (l — 5),x m ], where 



7T" 



2m 2 pq V — v 



i 

V+v 



Proof: Differentiate the explicit formula (JZJ) with respect to x to determine 
the critical points and sets Uq, U\, U%. The calculation of the resonance 
point in U\, U% requires to find two points in some neighborhood such that 
the derivative is strictly monotone on the interval between them, and has 
different signs at the extremities. □ 




Fig. 1: a) Typical form of the domains Uq, U\ and U^- Typical form of r] m (e) 
when {p,q) belongs Ui (b), U\ (c) and Uq (d). 



Remarks: 
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1. The optimal tuning rule can be rewritten in the form 

7T 



m(e) 



V V+v 

e 2 = . 



\J2pq V V — v 
The maximal value of amplification is found as 

2. We also see that the spectral power amplification as a measure of quality 
of stochastic resonance allows to distinguish a unique time scale, find its 
exponential rate (A = (V + v)/2) together with the pre-exponential factor. 
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d) 



Fig. 2: a) i] m (e) for p = q = 0.5, m = 500, v — 2, V — 4. Numerical 
simulations of X rn (k) for b) e = 0.65, c) £ = 0.9 and d) £ = 0.4. 
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